Modelling Current Filaments

In [1]:
import plotly.graph_objects as go
import plotly.express as px
import numpy as np
import random
from scipy.spatial.transform import Rotation as R

Define a filament object and an init function that creates one given a start point, azimuth, and phi - the angle from the +z axis

Code below used from previous work on PEM and XYZ Plate visualizer for calculating the BField

In [2]:
class Filament:
    def __init__(self, start, end):
        self.start = start
        self.end = end
    def __repr__(self):
        return "Start: {}, End: {}".format(self.start, self.end)
    
    def BField(self, array_x, array_y, array_z, current):
        """
        array_x, array_y, array_z must be ndarrays of equal shape with at least degree 3
        I is a current value in amperes
        """
        u0 = 1.25663706e-6
        loop_array = np.array([self.start, self.end])
        loop_diff = np.append(np.diff(loop_array, axis=0), [loop_array[0] - loop_array[-1]], axis=0)[:-1]
        loop_shift = np.append(loop_array[1:], [loop_array[0]], axis=0)[:-1]
        loop_array = loop_array[:-1]
        
        point = np.array([array_x, array_y, array_z])
        
        points = np.expand_dims(point, axis=1)
        #### This tile function is fixed to (1, 2, 1, 1, 1) because there are only 2 points in a filament, 
        #### not a closed loop or multi segment
        points = np.tile(points, (1, 2, 1, 1, 1))
        points = points.T
        
        # Vectors from points to endpoints
        AP = points - loop_array
        BP = points - loop_shift
        
        # Field Math
        r1 = np.sqrt((AP ** 2).sum(-1))[..., np.newaxis].T.squeeze().T
        r2 = np.sqrt((BP ** 2).sum(-1))[..., np.newaxis].T.squeeze().T
        Dot1 = np.multiply(AP, loop_diff).sum(-1)
        Dot2 = np.multiply(BP, loop_diff).sum(-1)
        cross = np.cross(loop_diff, AP)
        CrossSqrd = (np.sqrt((cross ** 2).sum(-1))[..., np.newaxis]).squeeze() ** 2
        top = (Dot1 / r1 - Dot2 / r2) * u0 * current
        bottom = (CrossSqrd * 4 * np.pi)
        factor = (top / bottom)
        factor = factor[..., np.newaxis]

        field = cross * factor
        field = np.sum(field, axis=3)
        field = field.T
        return field[0], field[1], field[2]
    
    def get_scatter3d(self, name, mode='lines+markers'):
        line = go.Scatter3d(x=[self.start[0], self.end[0]],
                            y=[self.start[1], self.end[1]],
                            z=[self.start[2], self.end[2]],
                            name=name,
                            mode=mode)
        return line
    
    def get_cones(self, name, x_array, y_array, z_array, current):
        U, V, W = U, V, W = newWire.BField(x_array, y_array, z_array, current)
        cones = go.Cone(x=x_array.flatten(), y=y_array.flatten(), z=z_array.flatten(), 
                        u=U.flatten(), v=V.flatten(), w=W.flatten(), 
                        sizemode='scaled',
                        sizeref=1,
                        anchor='tail',
                        name=name)
        return cones
In [3]:
# initialization function that takes in different parameters instead of a start, stop

def new_filament(start, length, deg_azimuth, deg_phi):
    """
    Create a new filament at P_start with length oriented in the 
    direction specified by azimuth, phi 
    """
    start = np.array(start)
    azi = deg_azimuth * np.pi / 180
    phi = deg_phi * np.pi / 180
    # Translate out of spherical coordinates
    x,y,z = np.sin(azi), np.cos(azi), np.cos(phi)
    dir_hat = np.array([x, y, z]) / np.linalg.norm(np.array([x, y, z]))
    direction = length * dir_hat
    endpoint = start + direction
    return Filament(start, endpoint)
In [4]:
# Testing method a
new_filament([0,0,0], length=1.732, deg_azimuth=45, deg_phi=45)
Out[4]:
Start: [0 0 0], End: [0.99997067 0.99997067 0.99997067]
In [5]:
# Testing method b
a = np.array([0,0,0])
b = np.array([1,1,1])
Filament(a, b)
Out[5]:
Start: [0 0 0], End: [1 1 1]

Draw a filament

In [6]:
a = np.array([0,0,0])
b = np.array([1,2,5])
newWire = Filament(a, b)
# Domain to draw the field over
X, Y, Z = np.meshgrid(np.linspace(-5,5,num=10),
                      np.linspace(-5,5,num=10),
                      np.linspace(-5,5,num=10))

line = newWire.get_scatter3d('NewLine')
cones = newWire.get_cones('Field', X, Y, Z, current=5)
fig = go.Figure(data=[cones, line], layout=dict(title=dict(text="B Field Plot")))
fig.show()
In [7]:
# Evaluate B Field at a point
i,j,k = 1,1,1

a,b,c = newWire.BField([[[i]]], [[[j]]], [[[k]]], current=5)
print(a.flatten(), b.flatten(), c.flatten())
[-1.14860606e-06] [1.53147475e-06] [-3.82868688e-07]

TESTING: Try summing some random current filaments and their corresponding fields

In [8]:
# Declarations
TX_CURRENT = 1000
NUM_LINES = 10
In [9]:
# Domain to draw the field over
x,y,z = np.meshgrid(np.linspace(-10,10,num=20),
                    np.linspace(-10,10,num=20),
                    np.linspace(-10,10,num=20))

filaments = []
for i in range(NUM_LINES):
    filaments.append(new_filament(start=[random.randint(-5,5),
                                         random.randint(-5,5),
                                         random.randint(-5,5)], 
                                  length=random.randint(1,10), 
                                  deg_azimuth=random.randint(0,359), 
                                  deg_phi=random.randint(0,180)))
    
fig = go.Figure(layout=dict(title=dict(text="Summed B Field Plot")))

Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for f in filaments:
    fx, fy, fz = f.BField(x, y, z, current=TX_CURRENT)
    fig.add_trace(f.get_scatter3d(str(random.randint(1,1000))))
    Bx += fx
    By += fy
    Bz += fz
cones = go.Cone(x=x.flatten(), y=y.flatten(), z=z.flatten(), 
                u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(), 
                sizemode='scaled',
                sizeref=1,
                anchor='tail',
                name='Compiled Field')
fig.add_trace(cones)
fig.show()
In [10]:
## Visualize with streamlines
figstream = go.Figure(layout=dict(title=dict(text="Summed B Field Plot - Streamline")))
for f in filaments:
    figstream.add_trace(f.get_scatter3d(str(random.randint(1,1000))))

# Initialize starting positions of streamlines
streamstarts = dict(x=[random.randint(-10,10) for i in range(100)],
                    y=[random.randint(-10,10) for i in range(100)],
                    z=[random.randint(-10,10) for i in range(100)])
streamlines = go.Streamtube(x=x.flatten(), y=y.flatten(), z=z.flatten(), 
                            u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(),
                            starts=streamstarts)
figstream.add_trace(streamlines)
figstream.show()

3D Brownian Motion Total Field of Resultant

Discretize $T$ into N intervals $\frac{T}{N-1}$ wide as $dt$

In [11]:
N = 100
T = 1.0
dt = T/(N-1)
"""
Create a line using a random walk algorithm.
dt = arbitary t Interval
N = number of lines to generate
"""
dX = np.sqrt(dt) * np.random.randn(1, N)
dY = np.sqrt(dt) * np.random.randn(1, N)
dZ = np.sqrt(dt) * np.random.randn(1, N)
X = np.cumsum(dX, axis=1)
Y = np.cumsum(dY, axis=1)
Z = np.cumsum(dZ, axis=1)

points = np.vstack((X, Y, Z)).T
In [12]:
f_list = []
for i, p in enumerate(points[1:]):
    f_list.append(Filament(points[i], p))

x,y,z = np.meshgrid(np.linspace(np.min(points[:,0]),np.max(points[:,0]),num=20),
                    np.linspace(np.min(points[:,1]),np.max(points[:,1]),num=20),
                    np.linspace(np.min(points[:,2]),np.max(points[:,2]),num=20))

fig = go.Figure(layout=dict(title=dict(text="Brownian Motion Simulation - BField")))
Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
for fil in f_list:
    fx, fy, fz = fil.BField(x, y, z, current=TX_CURRENT)
    fig.add_trace(fil.get_scatter3d(str(random.randint(1,1000)), mode='lines'))
    Bx += fx
    By += fy
    Bz += fz
    
cones = go.Cone(x=x.flatten(), y=y.flatten(), z=z.flatten(), 
                u=Bx.flatten(), v=By.flatten(), w=Bz.flatten(), 
                sizemode='scaled',
                sizeref=1,
                anchor='tail',
                name='Compiled Field')
fig.add_trace(cones)
fig.show()

This doesn't seem that useful in itself. But lets try simulating borehole paths between two electrodes and profile the Primary Field rotated to the hole geometry (AUV)

Define some survey geometry:

In [13]:
# Collar points of the boreholes 
collar1 = np.array([550,550,0])
collar2 = np.array([700,850,0])
# e1 and e2 are the real space representations of the electrodes
e1 = np.array([600,600,-100])
e2 = np.array([700,800,-150])

# Add an extra borehole that we will "read" stations down
reading_collar = np.array([650,650,100])
reading_bottom = np.array([650,740,-200])

Generate some stations so we can profile the primary down the hole...

In [14]:
NUM_STATIONS = 40
# Get some stations between these points
v = reading_bottom - reading_collar
print("Borehole length: {}".format(np.linalg.norm(v)))

# Generate some stations
v_hat = v / np.linalg.norm(v)
stations = np.array([reading_collar + v_hat * x for x in np.linspace(0, np.linalg.norm(v), num=NUM_STATIONS)])
station_interval = np.linalg.norm(stations[1]-stations[0])
print("Length intervals: {}".format(station_interval))
station_names = np.linspace(0, np.linalg.norm(v), num=NUM_STATIONS)
stations
Borehole length: 313.2091952673165
Length intervals: 8.031005006854258
Out[14]:
array([[ 650.        ,  650.        ,  100.        ],
       [ 650.        ,  652.30769231,   92.30769231],
       [ 650.        ,  654.61538462,   84.61538462],
       [ 650.        ,  656.92307692,   76.92307692],
       [ 650.        ,  659.23076923,   69.23076923],
       [ 650.        ,  661.53846154,   61.53846154],
       [ 650.        ,  663.84615385,   53.84615385],
       [ 650.        ,  666.15384615,   46.15384615],
       [ 650.        ,  668.46153846,   38.46153846],
       [ 650.        ,  670.76923077,   30.76923077],
       [ 650.        ,  673.07692308,   23.07692308],
       [ 650.        ,  675.38461538,   15.38461538],
       [ 650.        ,  677.69230769,    7.69230769],
       [ 650.        ,  680.        ,    0.        ],
       [ 650.        ,  682.30769231,   -7.69230769],
       [ 650.        ,  684.61538462,  -15.38461538],
       [ 650.        ,  686.92307692,  -23.07692308],
       [ 650.        ,  689.23076923,  -30.76923077],
       [ 650.        ,  691.53846154,  -38.46153846],
       [ 650.        ,  693.84615385,  -46.15384615],
       [ 650.        ,  696.15384615,  -53.84615385],
       [ 650.        ,  698.46153846,  -61.53846154],
       [ 650.        ,  700.76923077,  -69.23076923],
       [ 650.        ,  703.07692308,  -76.92307692],
       [ 650.        ,  705.38461538,  -84.61538462],
       [ 650.        ,  707.69230769,  -92.30769231],
       [ 650.        ,  710.        , -100.        ],
       [ 650.        ,  712.30769231, -107.69230769],
       [ 650.        ,  714.61538462, -115.38461538],
       [ 650.        ,  716.92307692, -123.07692308],
       [ 650.        ,  719.23076923, -130.76923077],
       [ 650.        ,  721.53846154, -138.46153846],
       [ 650.        ,  723.84615385, -146.15384615],
       [ 650.        ,  726.15384615, -153.84615385],
       [ 650.        ,  728.46153846, -161.53846154],
       [ 650.        ,  730.76923077, -169.23076923],
       [ 650.        ,  733.07692308, -176.92307692],
       [ 650.        ,  735.38461538, -184.61538462],
       [ 650.        ,  737.69230769, -192.30769231],
       [ 650.        ,  740.        , -200.        ]])

Set a hemispherical domain for each electrode, facing the other one to create a domain to iterate over for brownian movement

Pseudocode to generate a point from $P_1$ to $P_2$ where $P_1\neq P_2$, $P_{1z}\neq P_{2z}$ and $P_1$ and $P_2$ are position vectors in R3

We need to check the tolerance when walking to $P_2$. We will see if the generated point $P_n$ to $P_2$ is within tolerance. We will use a tolerance value where $|P_n - P_2| < Tolerance$ because generated floating point math is basically impossible to intersect exactly

Pseudocode:

function get_random_line:

create_a_ellipsoid_cone_domain

generate_a_point_in_domain

if point_in_tolerance (segment_length)
    just append last point
    return
else:
    append_generated_point
    recurse via get_random_line(generated_point)
In [15]:
# RANGE LIMITS
AZIMUTH_RANGE = 85
PHI_RANGE = 85

PHI_RANGE *= np.pi/180
AZIMUTH_RANGE *= np.pi/180

def vector_ranges(start_p, end_p):
# Vector from e1 to e2
    v = end_p - start_p
    azi_rad = np.arctan2(v[0], v[1])
    phi_rad= np.arccos(v[2] / np.linalg.norm(v))

    # Because arctan resolves to values +-pi, we have to 
    # adjust for the negatives to get proper azimuths
    if azi_rad < 0:
        azi_rad += 2 * np.pi

    azi_range = ((azi_rad + AZIMUTH_RANGE) % (2*np.pi),
                 (azi_rad - AZIMUTH_RANGE) % (2*np.pi))

    # We will bind phi to the domain [0, 180] 
    # so that it doesn't wrap around
    phi_range = (min(phi_rad + PHI_RANGE, np.pi),
                 max(phi_rad - PHI_RANGE, 0))
    return azi_range, phi_range

a_r, p_r = vector_ranges(e1, e2)

print('Azimuth Range: {}'.format(np.array(a_r) * (180 / np.pi)))
print('Phi Range: {}'.format(np.array(p_r) * (180 / np.pi)))
Azimuth Range: [111.56505118 301.56505118]
Phi Range: [180.          17.60438265]

Now we have the elliptic cone domain to call every step in the random walk, recall this function every step to regenerate the azimuth and phi ranges

In [16]:
def generate_random_points(start_p, target_p, segment_len):
    azi_range, phi_range = vector_ranges(start_p, target_p)
    
    azi_range = np.array(azi_range) * (180 / np.pi)
    
    # Catch the wraparound for azimuth across North
    if abs(azi_range[0] - azi_range[1]) > 180:
        azi_range = [max(azi_range) - 360, min(azi_range)]
        
    phi_range = np.array(phi_range) * (180 / np.pi)
    newLine = new_filament(start=list(start_p), 
                             length=segment_len, 
                             deg_azimuth=random.uniform(azi_range[0], azi_range[1]), 
                             deg_phi=random.uniform(phi_range[0], phi_range[1]))

    distance = np.linalg.norm(newLine.end - target_p)
    # Our tolerance distance will be the length of one segment
    if distance < segment_len:
        yield Filament(start_p, target_p)
    else:
        yield newLine
        yield from generate_random_points(newLine.end, target_p, segment_len)
In [17]:
FILAMENT_LENGTHS = 20
# Test by generating paths between the two electrodes, 
# with Filaments of length FILAMENT_LENGTHS
alistofwires = list(generate_random_points(e1, e2, FILAMENT_LENGTHS))
"Number of points in generated line: {}".format(len(alistofwires))
Out[17]:
'Number of points in generated line: 19'

Plotting the created line

In [18]:
fig = go.Figure(layout=dict(title=dict(text="Random motion between two points")))
fig.add_trace(go.Scatter3d(x=[e1[0], e2[0]],
                           y=[e1[1], e2[1]],
                           z=[e1[2], e2[2]],
                            name="Line Between Electrodes",
                            mode='lines+markers'))
for f in alistofwires:
    fig.add_trace(f.get_scatter3d(str(random.randint(1,1000)), mode='lines'))


fig.show()

We now have an algorithm to simulate possible current pathing. By iterating through random paths and generating theoretical response curves, we can generate loss functions which can be minimized to find a fit...

In [19]:
TX_CURRENT = 2
In [20]:
# Define a new function that sums the fields in a list of filaments 
def return_summed_BField(listofwires, x, y, z, current):
    Bx, By, Bz = np.zeros(x.shape), np.zeros(y.shape), np.zeros(z.shape)
    for fil in listofwires:
        fx, fy, fz = fil.BField(x, y, z, current=current)
        Bx += fx
        By += fy
        Bz += fz
    return np.array([Bx, By, Bz])
In [21]:
# Get station coordinates
stnx = np.array(stations[:,0])
stny = np.array(stations[:,1])
stnz = np.array(stations[:,2])
In [22]:
BField_stations = [return_summed_BField(alistofwires, 
                   np.array([[[sx]]]), 
                   np.array([[[sy]]]), 
                   np.array([[[sz]]]), 
                   TX_CURRENT) \
                   for sx, sy, sz in zip(stnx, stny, stnz)]
# Cleanup
for i in range(len(BField_stations)):
    BField_stations[i] = BField_stations[i].flatten()

BField_stations = np.array(BField_stations)
In [23]:
# Borehole vector
v
Out[23]:
array([   0,   90, -300])
In [24]:
# For our example, dip is constant - straight line borehole
BH_dip = np.arctan2(v[2], np.linalg.norm(v[:2])) * -180/np.pi
BH_dip
Out[24]:
73.30075576600638
In [25]:
# Borehole azimuth
BH_azimuth = np.arctan2(v[0], v[1]) * 180/np.pi
BH_azimuth
Out[25]:
0.0
In [26]:
print('XYZ Primary BField at Stations')
BField_stations
XYZ Primary BField at Stations
Out[26]:
array([[ 1.21676629e-09, -6.91857212e-10, -2.59779689e-10],
       [ 1.29476203e-09, -7.42793594e-10, -2.81180352e-10],
       [ 1.37981856e-09, -7.99606501e-10, -3.04989294e-10],
       [ 1.47270306e-09, -8.63219201e-10, -3.31495535e-10],
       [ 1.57426875e-09, -9.34733525e-10, -3.61008561e-10],
       [ 1.68546238e-09, -1.01546863e-09, -3.93850241e-10],
       [ 1.80733180e-09, -1.10700820e-09, -4.30340385e-10],
       [ 1.94103327e-09, -1.21125763e-09, -4.70772905e-10],
       [ 2.08783901e-09, -1.33051264e-09, -5.15378402e-10],
       [ 2.24914521e-09, -1.46754084e-09, -5.64267655e-10],
       [ 2.42648194e-09, -1.62567758e-09, -6.17348823e-10],
       [ 2.62152717e-09, -1.80893713e-09, -6.74209404e-10],
       [ 2.83612882e-09, -2.02214053e-09, -7.33952044e-10],
       [ 3.07234083e-09, -2.27106269e-09, -7.94971126e-10],
       [ 3.33248208e-09, -2.56260438e-09, -8.54653913e-10],
       [ 3.61922977e-09, -2.90500208e-09, -9.08983963e-10],
       [ 3.93576311e-09, -3.30810155e-09, -9.52010754e-10],
       [ 4.28597772e-09, -3.78373955e-09, -9.75118982e-10],
       [ 4.67480024e-09, -4.34629950e-09, -9.65967576e-10],
       [ 5.10864980e-09, -5.01351740e-09, -9.06845606e-10],
       [ 5.59612857e-09, -5.80757871e-09, -7.71966518e-10],
       [ 6.14909344e-09, -6.75638222e-09, -5.22822461e-10],
       [ 6.78439619e-09, -7.89434171e-09, -1.00040060e-10],
       [ 7.52682193e-09, -9.26071341e-09,  5.90866949e-10],
       [ 8.41408298e-09, -1.08898693e-08,  1.70350930e-09],
       [ 9.50449612e-09, -1.27795869e-08,  3.48647320e-09],
       [ 1.08832118e-08, -1.48086799e-08,  6.31775136e-09],
       [ 1.26364195e-08, -1.65744008e-08,  1.06660430e-08],
       [ 1.46717301e-08, -1.72437171e-08,  1.68048281e-08],
       [ 1.61557600e-08, -1.58470022e-08,  2.42706372e-08],
       [ 1.48736547e-08, -1.19894885e-08,  3.16700202e-08],
       [ 8.71705808e-09, -5.71646584e-09,  3.62569769e-08],
       [ 2.24438929e-10,  1.49671794e-09,  3.45636486e-08],
       [-5.57731298e-09,  6.46998597e-09,  2.80753193e-08],
       [-7.93474884e-09,  8.39294290e-09,  2.11904013e-08],
       [-8.39723671e-09,  8.46845933e-09,  1.58175329e-08],
       [-8.08704933e-09,  7.77705768e-09,  1.19741958e-08],
       [-7.50950897e-09,  6.86934397e-09,  9.24565533e-09],
       [-6.86434158e-09,  5.97331616e-09,  7.27218026e-09],
       [-6.22867637e-09,  5.16911639e-09,  5.81050059e-09]])

Now we rotate from XYZ to AUV

In [27]:
# Rotate the theoretical values into the same frame of reference used with boreholes
rotatedBField = R.from_euler('Z', -90, degrees=True).apply(BField_stations)

# Rotate the theoretical values into the hole coordinate system
auvField = R.from_euler('YZ', [90 - BH_dip, BH_azimuth], degrees=True).apply(rotatedBField)
print('Rotated AUV Primary BField at Stations')
auvField
Rotated AUV Primary BField at Stations
Out[27]:
array([[-7.37326167e-10, -1.21676629e-09, -5.00201072e-11],
       [-7.92263808e-10, -1.29476203e-09, -5.58817634e-11],
       [-8.53522153e-10, -1.37981856e-09, -6.23615251e-11],
       [-9.22068582e-10, -1.47270306e-09, -6.94709250e-11],
       [-9.99047387e-10, -1.57426875e-09, -7.71897871e-11],
       [-1.08581458e-09, -1.68546238e-09, -8.54473499e-11],
       [-1.18397895e-09, -1.80733180e-09, -9.40948665e-11],
       [-1.29545000e-09, -1.94103327e-09, -1.02866343e-10],
       [-1.42249287e-09, -2.08783901e-09, -1.11322988e-10],
       [-1.56779031e-09, -2.24914521e-09, -1.18775635e-10],
       [-1.73451060e-09, -2.42648194e-09, -1.24177914e-10],
       [-1.92638017e-09, -2.62152717e-09, -1.25981231e-10],
       [-2.14775892e-09, -2.83612882e-09, -1.21940756e-10],
       [-2.40371681e-09, -3.07234083e-09, -1.08859178e-10],
       [-2.70011283e-09, -3.33248208e-09, -8.22510332e-11],
       [-3.04368197e-09, -3.61922977e-09, -3.59025279e-11],
       [-3.44214490e-09, -3.93576311e-09,  3.87150631e-11],
       [-3.90436357e-09, -4.28597772e-09,  1.53254967e-10],
       [-4.44056864e-09, -4.67480024e-09,  3.23670837e-10],
       [-5.06265892e-09, -5.10864980e-09,  5.72023067e-10],
       [-5.78447449e-09, -5.59612857e-09,  9.29385642e-10],
       [-6.62167241e-09, -6.14909344e-09,  1.44065905e-09],
       [-7.59015429e-09, -6.78439619e-09,  2.17260140e-09],
       [-8.70037036e-09, -7.52682193e-09,  3.22699431e-09],
       [-9.94110323e-09, -8.41408298e-09,  4.76084689e-09],
       [-1.12387936e-08, -9.50449612e-09,  7.01162295e-09],
       [-1.23687504e-08, -1.08832118e-08,  1.03065512e-08],
       [-1.28105318e-08, -1.26364195e-08,  1.49788353e-08],
       [-1.16876537e-08, -1.46717301e-08,  2.10510518e-08],
       [-8.20455898e-09, -1.61557600e-08,  2.78006569e-08],
       [-2.38353388e-09, -1.48736547e-08,  3.37795320e-08],
       [ 4.94298441e-09, -8.71705808e-09,  3.63704999e-08],
       [ 1.13653871e-08, -2.24438929e-10,  3.26758924e-08],
       [ 1.42645063e-08,  5.57731298e-09,  2.50321420e-08],
       [ 1.41279983e-08,  7.93474884e-09,  1.78850289e-08],
       [ 1.26564476e-08,  8.39723671e-09,  1.27170549e-08],
       [ 1.08898301e-08,  8.08704933e-09,  9.23447838e-09],
       [ 9.23635773e-09,  7.50950897e-09,  6.88184024e-09],
       [ 7.81104485e-09,  6.86434158e-09,  5.24906563e-09],
       [ 6.62075061e-09,  6.22867637e-09,  4.08011553e-09]])
In [28]:
# Geometry Plot
x, y, z = [e1[0], e2[0]], [e1[1], e2[1]], [e1[2], e2[2]]
BH1 = [[e1[0], collar1[0]], [e1[1], collar1[1]], [e1[2], collar1[2]]]
BH2 = [[e2[0], collar2[0]], [e2[1], collar2[1]], [e2[2], collar2[2]]]
BHReading = [[reading_collar[0], reading_bottom[0]],
             [reading_collar[1], reading_bottom[1]],
             [reading_collar[2], reading_bottom[2]]]
fig = go.Figure()
fig.update_layout(title='Borehole MMR Profile Testing - Geometry')
fig.add_traces(go.Scatter3d(x=BH1[0], y=BH1[1], z=BH1[2], name='BH1'))
fig.add_traces(go.Scatter3d(x=BH2[0], y=BH2[1], z=BH2[2], name='BH2'))
fig.add_traces(go.Scatter3d(x=stations[:,0], 
                            y=stations[:,1], 
                            z=stations[:,2], 
                            name='Survey BH',
                            hovertext=station_names,
                            mode='lines'))
for f in alistofwires:
    fig.add_trace(f.get_scatter3d("Simulated Current Path", mode='lines'))
fig.show()

# AUV Profiling
Afig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"Z\"")))
Afig.add_trace(go.Scatter(x=station_names, y=auvField[:,0],
                         mode='lines+markers',
                         name='A Field'))
Afig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,0],
                         mode='lines+markers',
                         name='Z - Unrotated BField'))

Afig.show()

Ufig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"X\"")))
Ufig.add_trace(go.Scatter(x=station_names, y=auvField[:,1],
                         mode='lines+markers',
                         name='U Field'))
Ufig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,1],
                         mode='lines+markers',
                         name='X - Unrotated BField'))

Ufig.show()

Vfig = go.Figure(layout=dict(title=dict(text="Primary Field Profile \"Y\"")))
Vfig.add_trace(go.Scatter(x=station_names, y=auvField[:,2],
                         mode='lines+markers',
                         name='V Field'))
Vfig.add_trace(go.Scatter(x=station_names, y=BField_stations[:,2],
                         mode='lines+markers',
                         name='Y - Unrotated BField'))

Vfig.show()

This should serve as a framework for simulating direct current pathways moving forward...

In [ ]: